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HYPERVELOCITY IMPACT ANALYSIS BY THE 


METHOD OF CHARACTERISTICS 

By Richard Madden 
Langley Research Center 

SUMMARY 

A computer method has been developed for the analysis of the initial hydrodynamic 
phase of hypervelocity impact. The method is based on solving the inviscid, non heat- 
conducting, unsteady fluid dynamic equations in cylindrical coordinates by utilizing an 
approach based on the method of characteristics. The method has been used to generate 
an operational computer program applicable to the impact of an infinitely long right cir- 
cular cylinder striking a semi-infinite target of the same material. The projectile is 
assumed to be traveling parallel to its axis of symmetry, and the target is initially at 
rest. Results are presented for the pressure and velocity fields at various times fol- 
lowing the impact. The pressure fields are compared with the pressure fields calculated 
from an existing hydrodynamic analysis. 

INTRODUCTION 

The early stages of hypervelocity impact are characterized by shock waves in both 
the projectile and target materials with associated high pressures and internal energies. 
These extreme conditions justify the use of a hydrodynamic model for this portion of the 
impact. The hydrodynamic model has been customarily analyzed by using large-scale 
computer programs based on the particle in cell, or P.I.C., approach (ref. 1) or on a con- 
tinuous Eulerian extension of the P.I.C. approach (ref. 2). These approaches divide the 
flow field into cells and determine average values of the dependent variables (pressure, 
internal energy, etc.) over the cells. A disadvantage of the averaging procedure is that 
it causes numerical diffusion because discontinuities, such as shocks and rarefactions, 
are smeared over a number of cells. Thus, the items of prime concern are treated in a 
numerically undesirable fashion. These disadvantages have been overcome by utilizing a 
numerical solution based on the method of characteristics. The practical implementation 
of this method to a numerical solution is not, however, straightforward for problems 
involving more than two independent variables. 

The purpose of the present paper is to describe the development of a practical 
numerical solution of the hydrodynamic equations by using the method of characteristics 



and to discuss some of the techniques which were utilized to obtain an operational com- 
puter program. The method determines values of the dependent variables directly, at 
discrete points in the flow field, and the major discontinuities present in the flow field 
are located at each time step in the numerical calculation. This knowledge is used to 
develop a procedure which insures that linear interpolation relations and spatial finite 
difference relations are not applied across a discontinuity. For the discontinuities con- 
sidered, the numerical diffusion is essentially eliminated. 

In the present approach it is not possible to make a priori statements concerning 
the stability of the numerical procedure. However, the bicharacteristics have been 
selected such that they best represent the conoid of dependence of the point under consid- 
eration; thus, the probability of encountering a numerical instability is reduced. Although 
the stability of the procedure has not been verified mathematically, the numerical results 
give no indication of stability difficulties. The actual FORTRAN computer program uti- 
lized in the present approach is discussed in reference 3. 

SYMBOLS 


A coefficient in equation of state 

a coefficient in equation of state 

B coefficient in equation of state 

b coefficient in equation of state 

c speed of sound 

E internal energy 

E s sublimation energy 

E* coefficient in equation of state 

F function defining regions 

h time step 

l Q - cos i// 


2 



m Q = sin i// 


p pressure 

R radius of projectile 

r radial coordinate 

t time 

T elapsed time since impact 

U shock velocity 

u radial velocity 

u velocity normal to shock 

V initial velocity of projectile 

v axial velocity 

v velocity tangential to shock 

z axial coordinate 

a coefficient in equation of state 

/3 coefficient in equation of state 

£ coordinate tangential to shock 



0 angle between normal to characteristic surface and r-axis 

e = e - ^ 

X multiplier 
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P = V - 1 


£ coordinate normal to shock 

p mass density 

p* undisturbed mass density 

xp angle between normal to shock and r-axis 

Subscripts: 

i ith bicharacteristic 

o point in new time plane 

p particle 

s conditions ahead of shock 

GOVERNING EQUATIONS 

The initial stages of hypervelocity impact are characterized by very high pressures 
and very high internal energies. These extreme conditions allow the early stages to be 
analyzed by using a hydrodynamic model for the penetration process. (See, for example, 
ref. 2.) The present paper considers the initial deformation phase of a semi-infinite 
medium impacted by an infinitely long right circular cylinder. The cylinder and half- 
space are assumed to be of the same material. The problem as described is governed 
by the following classical, inviscid, non heat-conducting, Eulerian fluid dynamic equations 
with cylindrical symmetry: 

conservation of mass 


Dt 


+ P 


/ 9u 9v u' 
\8r 8z r 


= 0 


(la) 


conservation of radial momentum 


pgl + |E = o 

Dt 9r 


(lb) 
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conservation of axial momentum 


pDv + 9£ =0 (lc) 

H Dt 8z 

conservation of entropy along a particle path 

5R- C 2D2=0 (Id) 

Dt Dt 


where 


— = — + u — + v — 

Dt 8t 8r 8z 

p is density, u and v are radial and axial velocities, respectively, p is pressure, 
c is speed of sound, r and z are radial and axial coordinates, respectively, and 
t is time. 


Initial and Boundary Conditions 

For the problem under consideration, the initial conditions are that the infinitely 
long cylindrical projectile is moving along its axis of symmetry at a velocity V and the 
semi-infinite target is at rest (fig. 1). The thermodynamic state variables (pressure and 
internal energy) are assumed to be zero everywhere and the density is assumed to be 
undisturbed everywhere. 

The boundary conditions are that the radial velocity is zero on the axis of symmetry 
and the pressure is zero on the free surfaces of the projectile and target. The additional 
boundary conditions required to obtain a solution are applied on the low-pressure side of 
the shocks. In the present problem these conditions are that for both shocks the pressure 
and internal energy are zero, the density is undisturbed, and the velocity is equal to zero 
in the target and is equal to the impact velocity in the projectile. The boundary conditions 
are illustrated in figure 1(b) which shows a section taken through the axis of symmetry 
of figure 1(a). 


Equation of State 

The Tillotson equations of state (for either the compression region or the expansion 
region) from reference 4 have been used to describe the thermodynamic behavior of the 
materials. These equations may be written as follows: 

Ep + A p. + B p2 /p > p^ w ith 0 < E < E g ^ (2a) 


P = /a + 


E 


■+ 1 


E.,r? 
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and 



(p < p* with E > E s ) (2b) 


where 



fj. = rj - 1, and E s is the sublimation energy. 


GENERAL METHOD OF ANALYSIS 


Since in many problems it is desirable to obtain a more refined numerical repre- 
sentation of the flow field than may be obtained from the usual P.I.C. type of approach 
(refs. 1 and 2), it was decided to attack the problem by using an approach based on the 
method of characteristics in three independent variables. The hyperbolic nature of the 
equations makes them ideally suited for solution by the method of characteristics. This 
method allows discontinuities to be preserved and additionally permits the calculation of 
dependent variables at discrete points in the flow field; 

Characteristic Equations 

The characteristic equations are determined from the original partial differential 
equations (eqs. 1(a) to 1(d)) in appendix A. The final form of these equations, which 
replace equations (1), is taken from the appendix as follows: 

Slope equations along 
Bicharacteristic s 


— = u + c cos 9 
dt 


Particle paths 


dr 

dt 


u 


— = v + c sin 9 
dt 


dz 

dt 


v 


(3a) 


(3b) 


Compatibility equations along 
Bicharacteristics 


dp + pc du cos 9 + pc dv sin 9 = -pc 2 


— sin 2 0 
3r 


9u , 3v\ 

^ + 3 rj 


sin 9 cos 9 + — cos 2 0 + — 
3z r 


dt 

(3c) 
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Particle paths 


dE = -2- dp (3d) 

P 2 

where 9 is the angle between the normal to the bicharacteristic and the r-axis. These 
equations are highly complex since in addition to the total differentials there are also 
spatial partial derivatives. The solution of these equations must therefore rely on 
numerical techniques based on finite difference forms of the characteristic equations. 

Characteristic Equations in Numerical Form 

The original characteristic equations (3) have been approximated by using finite 
difference relations as follows: 

Slope equations along 
Bicharacteristics 


r j = r G - h^Uj + c^ cos 0^ 

(4a) 

Zj = z Q - h(vi + Ci sin 9 

(4b) 

Particle paths 


r p = r Q - hu p 

(4c) 

z p = z Q - hv p 

(4d) 

Compatibility equations along 
Bicharacteristics 


Po - Pi + Pi c i(u 0 " u i) cos % + Pi c i( v o ' v i) sin e i = -Pi c i 2h 

(If), ain2e ' 

• (£ + i?)i sln cos e i + (Hr), cos29 i + 

(4e) 

Particle paths 


E ° ‘ E P = 7^( p o “ p p) 

Pp 

(4f) 
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where the subscript i refers to values in the previous time plane associated with the 
bicharacteristic used and the subscript o refers to the point in the new time plane at 
which the solution is desired. The variables with subscript p refer to the particle 
located at ^rp,Zp^ in the previous time plane and which now resides at (r 0 ,z Q ). (See 
fig. 2.) Note that equations (4e) contain only three unknowns (p Q , u Q , and v 0 ) at some 
new point in time and that a numerical solution of equations (4e) along three bicharacter- 
istics is both necessary and sufficient for a unique numerical solution of these unknowns. 
These equations are solved for p Q , u Q , and v 0 to yield values on a spatial grid super- 
imposed on the flow field at each increment in time. 

Approximate Initial Conditions 

At the initial instant of impact (T = 0), conditions are specified on the surface; that 
is, z = 0, r > R. If the bicharacteristic slope and compatibility conditions could be inte- 
grated exactly, these initial data would be sufficient to determine the solution for any 
(r,z,t). The numerical procedure for the method of characteristics problem discussed 
herein utilizes more than two bicharacteristics and consequently requires initial data 
specified over some volume in the r,z space. 

Thus, some procedure is required to convert the initial conditions on the surface 
at T = 0 to a spatial set of initial conditions. This conversion was accomplished by 
advancing out the solution for some small time by an approximate technique. The 
resulting solution is then used as initial data for the characteristic solution in three inde- 
pendent variables. In this analysis the initial data are generated by considering the 
impact to be one -dimensional in the region directly beneath the projectile and in the com- 
pressed region of the projectile. These values are obtained directly from the solution of 
a one -dimensional impact. Values in the region r > R near the corner of the projectile 
are obtained by allowing the shock to advance spherically from the corner with a radius 
equal to the maximum depth of penetration of the target shock. The velocity of each par- 
ticle in this region is assumed to be equal to the particle velocity in the one -dimensional 
region; however, it is now directed radially from the point (R,0) to the particle. 

The rarefaction proceeding in from the free edges of the projectile is also included. 
The location of the rarefaction is determined from a spherical wave traveling at the 

speed of sound in the compressed material and whose center is at Values for 

all variables on the shocks, the shock location, and the rarefaction location are now known. 
Values for the interior points between the axis of symmetry and the rarefaction are the 
one-dimensional values, and points between the rarefaction and the free boundary are 
assigned values obtained by linear interpolation between the values on the free boundary 
and the one -dimensional values on the rarefaction. The values on the free boundary are 
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found by setting the pressure and internal energy equal to zero and the density equal to 
the undisturbed density p*. 

The present procedure for establishing initial conditions conserves neither 
momentum nor energy. However, if the radius of the projectile is large and the time 
step is small, the error introduced by the scheme is negligible. 

Approximate Boundary Conditions 

The boundary conditions are taken to be that the radial velocity is zero on the axis 
of symmetry and that the pressure vanishes along the original free surfaces of both the 
projectile and the target (r = R, z < 0 and z = 0, r ^ R). These pressure boundary 
conditions are not strictly correct because there is motion of the free surfaces and, con- 
sequently, the actual free surface does not coincide with the original free surface except 
at T = 0. This procedure was used for Convenience; the actual condition of a moving 
free surface could be handled by using the same location techniques as used for the shocks 
at the expense of a more complicated computer program. However, for the early stages 
of hypervelocity impact this approximation appears to be reasonable. 

NUMERICAL ANALYSIS 
Topological Aspects 

The numerical solution of equations (4) is complicated because physical disconti- 
nuities (shocks and rarefactions) appear in the flow field. Finite difference approxima- 
tions are meaningful only when the differenced function is continuous over the finite dif- 
ference interval. Therefore care must be taken to insure that finite difference relations 
are not applied across discontinuities. However, the discontinuities (shocks and rare- 
factions) divide the flow field into specific identifiable regions within which a numerical 
solution is valid. Thus all calculations which relate to a point in a given region may be 
required to use information from points within that region. Consequently, numerical 
diffusion (smearing out of discontinuities) is essentially eliminated for those disconti- 
nuities considered herein. The procedure for minimizing numerical diffusion is dis- 
cussed more generally in the following sections. 

Discontinuities and Particle Curves 

The procedure for minimizing diffusion requires locating the discontinuities at each 
time cycle. The discontinuities utilized in this analysis are as follows: 

(1) Projectile shock 

(2) Projectile -shock particle curve 
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(3) Target shock 

(4) Target -shock particle curve 

(5) Rarefaction 

(6) Rarefaction particle curve 

It is obvious that the shocks and rarefaction are locations of discontinuities. The addi- 
tional discontinuities which are defined herein as particle curves are the locus of all 
particles in the new time plane which were located on a given discontinuity in the old time 
plane. These curves must be defined since all particles which lie between this curve and 
the discontinuity in the new time plane have undergone discontinuous behavior during the 
time step. Consequently, the finite difference relations (4c), (4d), and (4f) written over 
the time step for any of these particles are meaningless. In this region, however, it is 
still possible to solve equations (4a), (4b), and (4e) to determine p 0 , u Q , and v 0 . 

This set of six discontinuities does not encompass all the discontinuities actually 
present in the physical situation. For instance, the reflection of a rarefaction from a 
shock or from the axis of symmetry is neglected. Discontinuities of this type have been 
neglected either because they are known to be weak or because they are not of major 
interest in hypervelocity impact problems. The neglect of these interactions as discon- 
tinuities does not remove them from the solution; however, it does allow them to be dif- 
fused throughout the flow field. 


Classification of Regions 

To categorize the behavior at a point in space, it is convenient to assign each point 
in the flow field to a particular region. The regions are shown in figure 3 and are defined 
as follows: 

Region I - Outside the target shock 
Region II - Outside the projectile shock 

Region in - Between the target shock and target-shock particle curve 
Region IV - Between the projectile shock and the projectile-shock particle curve 
Region V - Between the rarefaction and rarefaction particle curve 
Region VI - Outside the projectile and target 

Region VII - Between the rarefaction particle curve, the particle curves for the 
two shocks, and the free boundary 

Region VIII - Between the rarefaction, the particle curves for the two shocks, 
and the axis of symmetry 
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Separation of the flow field in this manner defines the regions in which values of the 
dependent variables may be calculated. 

In regions I, n, and VI, calculations are unnecessary since the material in these 
regions is undisturbed. In regions HI, IV, and V, calculations for density and internal 
energy may not be performed since the finite difference relation for the particle path 
equation is not applicable over the entire time step. Calculations for pressure and veloc- 
ity may be performed since the bicharacteristics are not discontinuous. However, for 
convenience, values of all dependent variables between the discontinuity and particle 
curves are obtained by interpolation between calculated values on the discontinuity and 
calculated values on an adjacent grid point behind the particle curve. In region VDI, cal- 
culations are unnecessary since the values of all dependent variables remain the same as 
the one -dimensional input values. In region VII, calculations are performed by using the 
procedures to be outlined in a subsequent section. 

Differentiation of Regions 

Each of the discontinuities at any time is specified by a given number of points. 
Between two successive points on the discontinuity a straight line may be constructed. 

This line has an equation of the form 

r = ajZ + 3-2 or z = a^r + a^ 

where aj . . . a^ are arbitrary constants. If these equations are written in a functional 
form F(r,z) such as 


F(r,z) = r - ajz - a 2 

(5a) 

F(r,z) = z - a0r - a^ 

(5b) 


certain properties can be recognized. In particular, on the discontinuity F is zero, 
and on one side of the discontinuity F is positive and on the other side F is negative. 

This property may best be illustrated by the following two examples. First, 
straight lines through rarefaction points are in the form of equation (5a). To determine 
whether a given point has been affected by the rarefaction, a straight line is constructed 
between the two closest rarefaction points and then F is evaluated. A value of F 
greater than zero indicates that the point has been affected by the rarefaction. Second, 
straight lines through shock points are in the form of equation (5b). To determine 
whether a given point lies outside the target shock, a straight line is constructed between 
the two closest target shock points and then F is determined. A value of F greater 
than zero indicates that the point is outside the target shock. 


11 



SOLUTION OF FINITE DIFFERENCE EQUATIONS 


The finite difference equations are solved at each time plane to yield values of the 
dependent variables on the projectile and target shocks and also at points on a rectangular 
grid superimposed on the flow field. The solution of the equations at various locations in 
the flow field requires the consideration of two classes of points: those points that lie 
sufficiently far away from a discontinuity (shock or rarefaction) and hence have a full 
conoid of dependence and those points that lie on or near a discontinuity and therefore 
have a partial conoid of dependence. The numerical representation of the equations and 
the procedures for handling these two classes of points are discussed in the following 
sections. A discussion of the stability of the numerical scheme is given in appendix B. 
The actual calculation procedure and a detailed analysis of the computer program with 
FORTRAN listing are presented in reference 3. 

Points Having a Full Conoid of Dependence 

For points having a full conoid of dependence, three bicharacteristics are required 
to approximate the conoid. The values for are chosen, equally spaced, to be 0, 

27r/3, and 47 t/ 3. The values of are chosen from stability considerations. (See 
appendix B.) The calculation requires the choice of a point (r 0 ,z 0 ) at which the solution 
is desired in the new time plane. The location (r^z^) in the previous time plane of the 
bicharacteristics which pass through (r 0 ,z Q ) may then be determined from equations (4a) 
and (4b). It is possible to find a point (rj,z^ for which u i? v^, and c^ satisfy equa- 
tions (4a) and (4b) since values of all dependent variables are specified on the grid and on 
the discontinuities in the previous time plane. However, the solution of equations (4a) 
and (4b) is complicated by the fact that there is no explicit relation tying (r^zjJ to 
( u i’ v i’ c i) * n a gi yen time plane. It is therefore necessary to choose some trial location 
for a bicharacteristic (rj,Zj) and use an iterative procedure (such as the Newton -Raph son 
method) to determine the actual coordinates rj and z ^ Once the coordinates r i 
and zj are established, interpolation is performed between the four surrounding grid 
points (or discontinuities) in the rz -plane to determine the values for all the required 
dependent variables at (ri,z^. Partial derivatives with respect to r and z at (r^,z^ 
are evaluated by using a forward or backward difference scheme. The three linear dif- 
ference equations (4e) are then solved simultaneously to give p Q , u Q , and v Q at 
(r 0 ,z 0 ). Equations (4c) and (4d) are then solved to determine the location in the previous 
time plane (rp,Zp^ of the particle which now resides at (r 0 ,z 0 ). Interpolation for depen- 
dent variables at ^rp,Zp^ is performed and the values for p Q and E Q are obtained from 
equation (4f) and the equations of state (2) by using an iterative procedure such as the 
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Newton -Raph son method. The speed of sound, at this point, may be obtained by differen- 
tiating the appropriate equation of state as follows: 


2 = 


9 P/e P 2 \ aE /p 


(6) 


Points Having a Partial Conoid of Dependence 

Points lying near a discontinuity .- Points having a partial conoid of dependence 
which lie near, but not on, a discontinuity are handled in almost the same manner as 
points having a full conoid of dependence. The only difference is that a new set of three 
bicharacteristics 0^ must be chosen to insure that none of the bicharacteristics cross 
the discontinuity. The bicharacteristics 6^ are established by a search technique which 
finds the intersection of the discontinuity and the conoid and then divides the remaining 
conoid into four equal angles. The computer program has a special routine which chooses 
the best set of bicharacteristics in this case. 

Points lying on a discontinuity .- The approach using three bicharacteristics is not 
required for points having a partial conoid of dependence which lie on a discontinuity 
since auxiliary relations or specified values of the dependent variables at these points 
supply sufficient information to form a complete set of equations. The particular pro- 
cedures applicable for each of these types of points are discussed separately. 

Points on the axis of symmetry: For points on the axis of symmetry, the radial 

velocity u vanishes (u = 0) and therefore only two bicharacteristics are required rather 
than the usual three. The two bicharacteristics 6^ chosen to be the best approximation 
of the conoid are 2n/3 and 4 tt/ 3. The solution of equations (4a) and (4b) for the two 
bicharacteristics is accomplished in the same manner as discussed for points having a 
full conoid of dependence. However, u G is set equal to zero in equations (4e) leaving 
only two variables p Q and v Q . The determination of p Q and E Q is performed in 
the same manner as explained in the section entitled "Points Having a Full Conoid of 
Dependence." 

Points on a free boundary: For points on a free boundary, the pressure p Q is 

equal to zero and again only two bicharacteristics are required. The choice of bicharac- 
teristics that lie in the flow field depends on whether the free boundary is along an r 
or z grid line. If r = R is the free boundary, the two bicharacteristics ^ and 

O 

5 7 r 

are chosen; whereas, if the boundary is along z = 0, the appropriate bicharacteristics 

O 

are 47r/3 and 57r/3. In the solution of equations (4e) the value of p Q is set equal to 
zero. To maintain p G = 0 on the free surfaces with the present equation of state, it was 
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assumed that the specific internal energy was zero (E = 0) and that density was undis- 
turbed (p = p*). 

Points on a shock wave: Points on a shock wave do not possess a full conoid of 
dependence since values of the dependent variables at these points can depend only on 
values at points encompassed by the shocks. The method employing three bicharacter- 
istics is not required inasmuch as additional relations are available from the conserva- 
tion laws expressed across the shocks. The procedure described for solving the equa- 
tions is similar to that developed in reference 5. 

The conservation relations across the shock and the equation of state provide four 
equations in the five unknowns and therefore only one bicharacteristic relation is required 
to form a complete set. 

For shock-wave calculations a set of local orthogonal coordinates £ is used at 
each point on the shock front. The coordinate £ is defined as the normal to the shock 
taken in the direction of motion, and the shock and the coordinate £ form a right-hand 
coordinate system with it. The £-axis forms an angle if/ with the original r-axis. The 
desired bicharacteristic is perpendicular to the shock and is 0 = 0 where 6 = 9 - ip. 
These relationships are shown in the following sketch: 



r 


Shock 


In this new set of axes, equation (3c) along a bicharacteristic becomes 


dp + pc cos 0 du + pc sin 9 dv = -pc^S dt 


( 7 ) 


where 


S = 9u sin 2£ _ 

n 


— + — Wn 9 cos 9 + — cos %9 + — 
a $} a? r 
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and 


0 = 0 -^ 


l 0 = COS I p 


u = l Q u + m Q v m 0 = sin \f/ 

v = Z 0 v + m Q u 


Note that the last term in S contains an unbarred u. The finite difference equivalent 
of equation (7) with 0 ~ \p (0=0) becomes 

Po "Pi + Pi c i( 5 o - 5 i) = -Pi c i 2h 

where the subscript 1 refers to the values at the point where the bicharacteristic 
( 8 = 0) intersects the previous time plane. The characteristic slope equations are most 
conveniently used in the form of equations (4a) and (4b) with 6 set equal to \p. The 
additional equations required are as follows: 

From conservation of mass and momentum across the shock, 



(“b - 5 s) 2 = ( p b _ p s)(^? ' p^) 

and 

u _ Pb% - Ps 5 s 


( 8 ) 


( 9 ) 


where the subscript s refers to conditions ahead of the shock and the subscript b 
refers to conditions just behind the shock. 

From the Rankine-Hugoniot relations, 

Eb-E s = i(Pb + P s )(^-i) (10) 

From the constancy of tangential velocity across the normal discontinuity (shock), 

v b = v s (11) 
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Equation (7a) is now combined with equation (8), after setting p g = 0, to obtain the fol- 
lowing equation which is more amenable to solution: 



If 1 1 \ 

2. 

/9v\ U 1 

P b + Pl c l 


= -Pi h 

( 8 ?) 1 + r l 


+ Pi + Pi c i^i 


( 12 ) 


where the positive square root is used in target shock calculations (u s = 0) and the neg- 
ative square root is used in projectile shock calculations (u s = v). 

Equations (10) and (12) are solved with p s = 0 by using an iterative procedure, 
such as the Newton -Raph son method, to obtain and E^. The pressure p^ in these 
equations is obtained from the equation of state (2a). The speed of sound may now be 
obtained from equation (6), and the velocity normal to the shock u^ is obtained from 
equation (8). 


RESULTS AND DISCUSSION 

The particular example chosen to illustrate the results of the characteristic 
approach is the impact of an infinitely long cylindrical aluminum projectile 2.5 cm in 
diameter on an aluminum half-space. The initial velocity of the projectile is taken to 
be 7.6 km/s and the half-space is assumed to be at rest. Initially, a 0.12 5 -cm -square 
grid is used; however, a grid change is made at T = 1.20 ps and thereafter the grid 
becomes a 0.25-cm-square grid. The Tillotson equations of state (ref. 4) are used to 
define the thermodynamic behavior of the aluminum. 

Pressure fields and velocity fields from the characteristic method are illustrated 
at T = 0.74 ps, T = 1.26 ps, and T = 1.50 ps. Also presented are a plot of pressure 
on the axis of symmetry as a function of time and a plot of discontinuity location as a 
function of time. The pressure plots are compared with the results from the well-known 
OIL computer program of reference 2. Results from the OIL program have been com- 
pared extensively with exact solutions and therefore are a logical choice for comparison 
herein. 


Pressure Fields 

Figure 4 gives the pressures on the axis of symmetry at a time before reflection of 
the rarefaction from the axis of symmetry (T = 0.74 ps) , obtained by both the method of 
characteristics and the OIL approach. Since the discontinuities have been tracked, the 
characteristic method predicts the exact hydrodynamic behavior, that is, a one- 
dimensional region bounded by two sharp shocks. The OIL approach, on the other hand, 
has numerical diffusion near the shocks. It is noted also that the OIL method tries to 
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compensate for the smearing by admitting a pressure directly behind the shock which is 
higher than the one-dimensional value. 

Figure 5 shows the pressure distributions in the entire flow field for both solutions 
at T = 0.74 /is. These distributions are compared by plotting isobars for p = 1.08, 

0.8, 0.5, and 0.2 Mb (108, 80, 50, and 20 GN/m2). (Note that 1 Mb = 100 GN/m2.) T he 
method of characteristics distributions are shown as solid lines and the OIL distributions 
as dashed lines. The locations of the two shocks and the rarefaction are clearly defined 
for the characteristic method. One difference immediately apparent between the two sets 
of isobars is that the isobars for the OIL approach begin and end on the axis of symmetry, 
whereas the isobars for the characteristic method begin and end on the shocks. This 
difference, of course, comes from the smearing of shocks in the OIL solution. The OIL 
approach additionally predicts two regions of pressures higher than the one -dimensional 
region bounded by the two shocks and the rarefaction (p = 1.08 Mb). The reduction of 
pressure in the rarefied region is predicted, similarly, by both methods. Agreement of 
the isobars for p = 0.8 Mb is good except near the shocks. The isobars for p = 0.5 
and 0.2 Mb also compare quite favorably with major differences noted near the free sur- 
face, where pressures are low. In general, the results of the two methods are in agree- 
ment. However, the absence of diffusion in the characteristic method results in a more 
accurate solution at any point. 

The pressure fields for both the characteristic and OIL approaches at T = 1.26 jus 
are shown in figure 6. At this time a portion of the rarefaction has been reflected from 
the axis of symmetry. Again, the trends for both approaches are in agreement. The 
difference between the isobars for p = 0.5 Mb in the central portion of the flow field is 
caused, in part, by the low pressure between the shocks computed by the OIL approach 
(fig. 4). 

Figure 7 illustrates the isobars as predicted by both approaches at T = 1.5 jus. 

The isobars again compare favorably. Considering the characteristic method results, the 
higher pressure isobars for the earlier flow stage were seen to follow the shape of the 
rarefaction. However, for this time and for succeeding stages in the flow the higher 
pressure isobars follow a different pattern. In particular, these isobars turn outward in 
the central region of the flow field; this indicates that the pressure is higher than might 
be expected from an extrapolation of the behavior of the isobars at an earlier stage in the 
flow. This higher pressure results from the fact that as the projectile shock progresses 
it is attenuated and, therefore, particles passing through it are not decelerated as much. 
These particles enter the area between the projectile and target shocks with high velocities 
and must again be decelerated in the flow field. This additional deceleration gives rise 
to the increase in pressure. 
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Figure 8 presents pressure profiles on the axis of symmetry at three different 
times. The profile at T = 1.26 ps illustrates a completely one -dimensional behavior 
and is similar to the profiles that would be obtained for times before reflection of the 
rarefaction from the axis of symmetry. The next type of profile is shown at T = 1.48 ps. 
At this time reflection of a portion of the rarefaction has occurred and gives an area in 
which the pressures are less than one -dimensional. However, the rarefaction has not 
completely reflected and therefore one -dimensional regions still remain adjacent to the 
shocks. The slope of the rarefaction near the projectile shock is greater than that near 
the target shock; this indicates that the rarefaction is much stronger near the projectile 
shock than near the target shock. This difference in strength is caused by the relative 
distances from the free surface. The profile at T = 1.75 ps illustrates the behavior 
following complete reflection of the rarefaction. At this time the entire profile is simi- 
lar to the center portion of the profile for T = 1.48 ps and additionally the shocks have 
been attenuated to a pressure below the one -dimensional value. 

Velocity Fields 

The velocity vector field is shown in figure 9 for T = 0.74 ps. Each arrow indi- 
cates both the magnitude and the direction of the velocity present at the tail of the vector. 
A vector of 7.6 km/s has been plotted below the flow field to establish the scale for the 
plot. The two shocks and the rarefaction are shown as solid lines. The vectors indicate 
that the area of the projectile unaffected by the shock has the original projectile velocity 
(7.6 km/s) and that the area between the shocks, rarefaction, and axis of symmetry is 
one -dimensional with a velocity equal to one -half the impact velocity. The velocity vec- 
tors are seen to turn radially outwards in the region outside the rarefaction. This 
change in direction results from the radial component imparted to the velocity vector 
when a point passes out of the rarefaction and from spherical divergence. The vectors 
indicate that there is a vortex-type flow as material slides along the target shock and 
moves upward through the original free surface, as evidenced by the negatively directed 
vectors near the target free surface. Some of this material eventually forms the lip 
which is commonly seen in impacted semi-infinite targets. The magnitude of the vectors 
illustrates that in the area near the rarefacted target shock, the velocity is less than the 
one -dimensional velocity and decreases radially along the shock. However, in the area 
of the projectile shock, velocities are above the one -dimensional value. This velocity 
increase results from the attenuation of the projectile shock strength and particles 
passing through the shock are not slowed greatly. In fact, near the projectile radius, 
where the projectile shock is strongly attenuated, some of the velocity vectors are 50 to 
60 percent greater than the one -dimensional value. 
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Figures 10 and 11 illustrate the velocity vector fields at T = 1.26 jtx s and 
T = 1.5 jus, respectively. Since the problem under consideration is the impact of an 
infinitely long projectile, the flow field is continually "fed" through the projectile shock 
and, consequently, the velocity vector plots have similar characteristics for all times. 
The only differences are in the continued "vortex" activity and the increased magnitude 
of the vectors near the projectile shock. In fact, the velocities near the intersection of 
the projectile shock and the free surface are just slightly reduced from the original 
projectile velocity. 


Discontinuity Locations 


In figure 12 the location of the discontinuities is shown as a function of time. The 
behavior of the rarefaction may be determined independently from that of the shocks 
since the rarefaction always proceeds toward the axis of symmetry as a toroidal segment 


with its center at the point 


¥) 


The shocks are seen to be straight before arrival of the rarefaction. After arrival 
of the rarefaction the strength of the target shock decreases and the shock normal turns 
radially. After the passage of less than 1 ps the shock velocity has been reduced to 
approximately the speed of sound in the uncompressed material. Since the velocity of the 
material in front of the target shock is zero, this shock always advances in the same 
direction. The projectile shock, however, behaves a little differently since the position 
of a projectile shock point is governed by both its local shock velocity and the impact 
velocity. At this impact velocity a projectile shock point, therefore, advances in the 
negative z-direction until its local shock velocity is attenuated to a value below the orig- 
inal projectile velocity. Thereafter, the point travels in a positive z-direction. This 
motion can be seen by observing the projectile shock point closest to the free surface. 


CONCLUDING REMARKS 


A numerical method of analysis for determining the initial hydrodynamic phase of 
the axisymmetric impact of an infinitely long right circular cylindrical projectile on a 
semi-infinite target has been developed and programed. This method, based on the 
method of characteristics in three independent variables, differs from other approaches 
currently in use, such as the particle in cell (P.I.C.) method, in that specific values may 
be obtained at any point in the flow field directly rattier than by integrating over some 
finite volume. The method keeps track of the locations of the major discontinuities, such 
as the shocks and rarefaction, at each time step and therefore essentially eliminates the 
numerical diffusion present in other numerical solutions to hydrodynamic behavior. Thus, 
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the shock-wave phenomenon, which is of prime interest, is not smeared out in the numer- 
ical calculation procedure. 

The pressure fields obtained by using this numerical method of analysis have been 
compared with pressure fields obtained from the computer program based on a continu- 
ous Eulerian extension of the particle in cell approach (OIL) and good agreement is indi- 
cated when smearing of discontinuities inherent in the OIL program is disregarded. Of 
the two programs, OIL takes approximately one-tenth of the computer time of the charac- 
teristic solution; the additional time required by the method of characteristics is justified 
by the increased point-to-point accuracy. The techniques used in developing the charac- 
teristic solution should be valuable in other areas of fluid dynamics as well as hyper - 
velocity impact. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., July 10, 1968, 

124-09-15-04-23. 
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APPENDIX A 


METHOD OF CHARACTERISTICS 

A brief description of the development of the characteristic slope and compatibility 
equations is presented. More general discussion of the method of characteristics and 
its applications to other fields is given in references 6 to 9. 

Characteristic Slope Equations 

The fluid dynamic equations (1) may be condensed to a convenient form by using 
"index notation" as follows: 


^Hjk 



+ bi 


= 0 


(Al) 


where <p j represent the dependent variables, x k represent the independent variables, 
and ajjk = (^j ) ^i = ^i(^j’ x k)' The characteristic slope equations for equa- 

tions (Al), most easily developed by changing the independent variables from x k to 
some arbitrary coordinate system /3j_, 02 , and 03 > are 


a ijk 


d _Vj 9/3 m 
90 m 9x k 


= 0 


(A2) 


With the assumption that values of all the dependent variables and their derivatives with 
respect to ^ an d 03 are specified on a surface 0^ = Constant, these transformed 
partial differential equations would be expected to yield the derivatives with respect to 
0j, if they exist. These derivatives are 


9(pj 90j 
&ijk 8^ 


-b. 


" ^jk 


9 0 a 


{a = 2,3) (A3) 


In developing the characteristic slope equations, however, it is desired to determine 
the conditions under which the derivatives normal to 0j do not exist (normal derivatives 
to the surfaces, 0j = Constant, are discontinuous). These surfaces are then called char- 
acteristic surfaces. The requirement for discontinuities in the derivatives with respect 
to 0j is then the vanishing of the determinant of the coefficients of the derivatives with 
respect to 0j in equations (A3); that is, 
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This condition always yields a set of slope equations; however, the /3j surfaces corre- 
spond to physical surfaces only when these slopes are real, that is, when the original 
equations are hyperbolic or parabolic. For equations (1), equation (A4) takes the form 


Dt 


0 

D/3! 

301 

3/3! 

n i. 

Dt 

P 9r 

9z 

3% 

0 

D^! 

0 

“IF 


p "dF 


3^ 

0 

0 

D/3! 

9z 



p "dT 

D/3! 

2 D % 
-c z — — 

0 

0 


Dt 


= 0 


(A5) 


Expanding the determinant gives 
H \ Dt , 

where 



+ 


9^ 

“IF 


8% 

9z 


= 0 


D/3. 8/3-. 80j "p. 

= — - + u + V 

Dt 9t 9r 9z 

Equation (A6) has three possible solutions: the trivial solution p = 0 and two nontrivial 
solutions 

^l)\o 

Dt / 


9 / 3 -, 


(A6) 


(A7) 


and 


D/3 


Dt 


11 +C 2 


9%\ 2 / 9%' 2 

9r / + \ 9z 


= 0 


(A8) 


Equation (A7) may be investigated as a linear equation; however, for the sake of consis- 
tency in the mathematical analysis it is treated as a nonlinear equation. By expressing 
equation (A7) as 


F = 


9% 9/3. 9)3., 

— - + u — -+v — -I =0 
9t 9r 9z 


(A9) 
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and applying the techniques from nonlinear partial differential equations (ref. 6), the 
characteristic slope equations become 

St- *(* ♦■*♦’*)" 


dz 

dv 


— = 2(0 t + u/3 r + v0 z )v 


dt_ 

dv 




+ V/3 Z ) 


where subscripts r, z, and t denote derivatives of and v is a parameter. 
Eliminating the parameter y gives 


dr 

dt 


u 



(A10) 


The characteristic slopes for equation (A7) are thus the particle paths which are surfaces 
of possible discontinuity in the entropy derivative. 

Equation (A8) may be analyzed in an analogous fashion by redefining F to be 

= 0 (All) 

The characteristic slope equations therefore are 



i? = ^:' 2 (' s t + u ' 5 r tv < s z ) u - 2c2,3 r 


dz 

du 


8j- = 2(/3 t + u£ r + v£ z )v -2 c 2 /3 z 
* z 


Eliminating the parameter v gives 

dr _ u c ^^r 

dt ” P t + u/3 r + v/3 z 


dz 

dt 


v 


% + u ^r + v ^z 
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Using equation (All) to reduce these equations yields 


d£=u±c 

dt 


{ 0 ? + / 3 Z 2 


dz 

dt 




(A12) 


= v ± c 




Equations (A12) may be reduced to a more tractable form by noting that the direction 
cosines of the projection on the rz-plane of the normal to the surface (3^ = Constant are 


cos 


M = 

[h’ z ) = 


•Jh 2 * 


\fe? * 


By defining 


the slope equations (A12) become 


-cos^j8j,r^ = cos 0 
-cos^jjZ^ = sin 9 


dr 

— = u ± c cos 9 
dt 


— = v ± c sin 9 


(A13) 




(A14) 


dt j 

In equations (A14) only the positive sign need be considered since the negative sign may 
be obtained by changing the reference for 9 by v. The characteristic slope equations 
are then 


dr 

— = u + c cos 9 
dt 


d 7 

— - = v + c sin 9 


(A15) 


dt 
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In equations (A15) a given value of the parameter 9 (for 0 S 8 S 2tt) defines one 
of the characteristic directions at a point. These characteristics are termed 
"bicharacteristics." Considering the entire range of 9 (for 0 S 9 = 2n), equations (A15) 
describe a general cone in space called the characteristic conoid (fig. 13). The family of 
bicharacteristics are the generators of the conoid. In the (t Q - h) plane, the conoid 
envelops an area termed "the domain of dependence" which encompasses all points which 
affect ^r 0 ,z 0 ) at t Q . The point (r 0 ,z O’*o) is dependent only on the conditions which exist 
inside and on the conoid of dependence. 

Compatibility Equations 

A compatibility equation is a partial or total differential equation in which none of 
the indeterminable derivatives normal to /3j exist — that is, it is an equation which is 
valid along the directions indicated by the characteristic slope equations. 

Equations (A10) are recognized to be particle path equations and therefore there 
are many possible compatibility equations for this slope, for example, equations (lb), 

(lc), and (Id). However, for the present problem the first law of thermodynamics taken 
along the isentropic particle path 


T ds = 0 = dE - dp (A 16) 

P l 

is the most useful (T being temperature and s being entropy). 

The compatibility equation corresponding to the bicharacteristics given by equa- 
tions (A15) is obtained by combining the transformed equations (1) in a manner such that 
the indeterminable derivatives (with respect to Pi) do not appear. This simplification is 
done most easily by multiplying equations (la), (lb), (lc), and (Id) by weighting factors 
A i, Ag, Ag, and A^, respectively, and summing. Relations between the A T s are found 
by equating to zero the coefficients of the derivatives with respect to j3j in the trans- 
formed equations. The derivatives can be written in the form 
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Applying this procedure gives 


D/3j 8j8, 8/3i 

Xa + Xo h X q — — 0 

4 Dt 1 9r 3 8z 

9/3i D/3, 

1 8r 2 Dt 

8j8i D/3i 

X, — — + Xo rrr - - 0 

1 dz <) Dt 


Xl ^i-X 4 c2^1=0 

1 Dt 4 Dt 

Solving for X 2 , X 3 , and X 4 in terms of Xj yields 


9% 

X - X 8r 

X 2- X l D/3l 

“Dt” 


X 3 - 


- X 


9% 

d Z 

1 

dT 


From equations (All) and (A13) 


x - Xl 

X 4--2 


9% 

8 r 

D/3j 

DT 


cos 6 


9% 

9z sin 0 
Dgl " c 
"dT 


(A17) 


(A18) 


(A19) 
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The negative sign results from the choice of reference for 9; therefore, 

Ai 

An = — cos 9 
* c 

*1 

An = — sin 9 
o c 


(A20) 


* 4 C 2 


By using equations (A20) for the A's, the compatibility equation is determined from 
the sum of the weighted equations (la) to (Id) as follows: 


5£+ p 

te + 

*l + *\ 


iElcos * ♦ il( 

pDv + 

9p> 

(sin 9 + - 


Dt P| 


9z r J 

c Dt 

9r J c \ 

' Dt 

9Z ; 

/ c2\Dt 

Dt / 


(A21) 

Along a bicharacteristic the time derivative becomes 


d _ 9 8 dr | 9 dz 

dt dt 9r dt dz dt 


= — + (u + c cos 9)-^-+ (v + c sin 9) 77 - 
9 t dr 9z 


(A22) 


Substituting equation (A22) into equation (A21) and reducing gives the final simplified form 
of the compatibility relation as 


dp + pc du cos 9 + pc dv sin 9 = -pc^ ^|H.j s in20 - + |^-jsin ® 


cos 9 + 



(A23) 
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NUMERICAL STABILITY 

The stability of numerical solutions to partial differential equations has been of 
much concern over the past few years. For hyperbolic equations, stability is usually 
based on the Von Neumann criterion (ref. 10) or the Courant-Friedrichs-Lewy (CFL) 
condition (ref. 11). The Von Neumann criterion, although very powerful, is not directly 
applicable since it relies on matrices and linear equations in one independent variable. 

The CFL condition is also difficult to apply since it requires that the domain of dependence 
of the difference equations encompass the domain of dependence of the original partial dif- 
ferential equations (fig. 14). It is possible, however, to use certain qualitative measures 
of stability based on the CFL philosophy. A brief discussion is presented of the approach 
used in the present study. 

Consider the domain of dependence of the partial differential equations in the 
(t Q - h) plane as illustrated in figure 15(a). If the conoid is now approximated by the 
three bicharacteristics indicated by the triangular symbols, a domain of dependence of 
the difference equations may be defined (ref. 10). The domain of dependence of the dif- 
ference equations is described by the convex hull of the base points used in the calculation. 
When only points 1, 2, and 3 are used, the convex hull is a triangle connecting the three 
points (dashed line). 

For stability, according to the CFL conditions the triangle must encompass the 
circle. The maximum probability of encompassing the circle will naturally occur when 
points 1, 2, and 3 are equally spaced around the circle; this indicates that the three 0^'s 
should be chosen at 2n/3 angular intervals. However, when only points 1, 2, and 3 are 
used, another problem arises from convergence considerations. Since convergence dic- 
tates that the points 1, 2, and 3 lie relatively close to the circle, the triangle will not 
encompass the circle and instabilities will exist. This problem must be eliminated by 
increasing the domain of dependence of the difference scheme by utilizing more points in 
the previous time plane. 

In the method utilized in the present paper, the points 1, 2, and 3 are obtained by 
interpolation in the calculation grid in the (t G - h) plane. The original three points are 
replaced by the 12 grid points (circular symbols) which surround points 1, 2, and 3 and 
thus the domain of dependence is increased to that shown in figure 15(b). This enlarged 
domain then gives a greater probability of stability. 

Since there has been no evidence of instability in the present results, the spreading 
of the bicharacteristics and the interpolation in the (t Q - h) plane appear to have assisted 
greatly in assuring numerical stability. 
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Figure 9.- Velocity vector field as predicted by the characteristic method at T = 0.74 ms for the example impact problem. 













OJ l> 



National Aeronautics and Space Administration 
Washington, D. C. 20546 


OFFICIAL BUSINESS 


I 1 4 '< 


i • . I 


FIRST CLASS MAIL 


if . \ ; ;v/ a i- . i / 
r , : M' ' 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS M 
SPACE ADMINISTRATION - 


/ ! i t 


dhctu a cTru . If Undeliverable ( Section 13 
POSTMASTER . p ostal ManuaI } Do Not R e P 


ft The aeronautical and space activities of the United States shall be 4 

conducted so as to contribute ... to the expansion of human knowl- 4 

edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof ” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


\ 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



